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ABSTRACT 

The star HD 69830 exhibits radial velocity variations attributed to three planets as well as 
infrared emission at 8 - 35/im attributed to a warm debris disk. Previo us studies have de- 
veloped models for the planet migrati on and ma ss growth dAlibert et al.ir2005l l2006) and the 
replenishment of warm grains ( Wvatt et al.ll20()71 In this paper, we perform n-body integra- 
f~| [ tions in order to explore the implications of these models for: 1) the excitation of planetary 

pi !■ eccentricity, 2) the accretion and clearing of a putative planetesimal disk, 3) the distribution of 

planetesimal orbits following migration, and 4) the implications for the origin of the infrared 
emission from the HD 69830 system. 

We find that: i) It is not possible to explain the observed planetary eccentricities (~ 0.1) 
^ I purely as the result of planetary perturbations during migration unless the planetary system is 

nearly face-on. However, the presence of gas damping in the system only serves to exacerbate 
the problem again, ii) The rate of accretion of planetesimals onto planets in our n-body sim- 
^ , ulations is significantly different to that assumed in the semi-analytic models, with our inner 

. planet accreting at a rate an order of magnitude greater than the outer ones, su ggesting that 

' one cannot successfully treat planetesimal accretion in the simplified manner of lAlibert et al.l 

, ([2006). iii ) We find that the eccentricity damping of planetesimals does not act as an insur- 

■ mountable obstacle to the existence of an excited eccentric disk: All simulations result in a 

significant fraction (~15%) of the total planetesimal disk mass, corresponding to ~ 25Me, 
remaining bound in the region ~l-9 AU, even after all three planets have migrated through 
the region iv) This swarm of planetesimals has orbital distributions that are size-sorted by 
I gas drag, with the largest planetesimals (~ l,000fem), which may contain a large proportion 

of the system mass, preferentially occupying the highest eccentricity (and thus longest-lived) 
orbits. Although such planetesimals would be expected to collide and produce a disk of warm 
k> \ dust, further work will be required to understand whether these eccentricity distributions are 

j_j ■ high enough to explain the level of dust emission observed despite mass loss via steady state 

d ' colUsional evolution. 

1 INTRODUCTION indicative of warm dust (other systems were observable only at 

70/jm, indicative of much cooler (T < lOOK) dust). 
Radial velocity surveys have discovered over 300 extrasolar plan- All three known planets orbit at less than 1 AU, so the putative 
ets and roughly 30 multiple planet systems (Butler et al. 2006; dust disk is exterior to the known planets, but still close enough 
http://www.exoplanets.org, http://www.exoplanet.eu). The multiple that gravitational perturbations from the planets might be effective 
planet systems are particularly valu able for gaining insi g hts into the at stirring a disk of planetesimals, resulting in collisions that could 
proces ses of planet formation (e.g.. lLee & Peai3 l l2002h : lFord eTal] generate dust that would give rise to the observed IR excess. Thus 
( |2005h ). Among the k nown multiple plan et systems, the Spitzer the HD 69830 system is particularly interesting due to a unique 
Legacy Program FEPS i Mever et al.ll2006l) has identified four sys- combination of radial velocity and IR excess observations that have 
terns that emit an excess of infrared radiation (relative to that ex- revealed multiple planets and a close-in warm dust disk, 
pected from the stellar photosphere), most likely due to reradiation 
of starlight absorbed by a dust disk. In these systems, the proper- 
ties of the dust disk can provide additional con straints on the sys- 1.1 Planets: Radial Velocity Observations and Theoretical 
tem's formation and dynamical evolution (e.g.. lMoro-Martm et al.l Models 

Based on HARPS radial-velocity measurements of HD 69830, 
One of th ese systems is HP 698 30, a bright KOV star located iLovis et al. 1 ( l2006h identified planets (Msini = 10.2, 11.8 & 



12.6 pc away jPerrvman et al. 1997). In addition to possessing 3 18.1Me respectively) orbiting close-in to the central star (semi- 



Neptune-mass planets jLovis et al.ll2006l) . this ex ceptional system major axes of 0.08, 0.19 & 0.63 AU). The planets' e ccentricities 
was th e only system out of 84 FKG stars studied bv lBeichman et al.l (0.10 ± 0.04, 0. 1 3 ± 0.06 & 0.07 ± 0.07 respectively, iLovis et al.) 
l l2005h which was found to display an excess of emission at 24yum, | |2006|) ) are large by solar system standards, but modest when com- 
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pared to those of extra-solar planets i lFord & Rasiol2007h . The host 
star is somewhat cooler {T^ff ^ 5385 K) and less massive than the 
Sun (M-i, ^ O.86M0). Otherwise, it is quite similar to the Sun, with 
a nearly solar met allicity ([Fe/H]= -0.05±0.02) and age (4-10 Gyr; 
iLovis et al. Since the known planets produce modest radial 

velocity perturbations (2.2 - 3.5m s"') just large enough for detec- 
tion, similar planetary systems could have often eluded detection by 
broad radial velocity planet searches on account of an insufficient 
number of velocity observations and/or velocity precision due to 
photon noise and/or radial velocity "jitter". Recent results from the 
HARPS search for southern extra-solar planets and its discovery of 
a planetary system with 3 Super-Earths ( IMavor et al. 2008) bolsters 
the view that similar planetary systems might be quite common. 

These radial velocity observations alone reveal an interesting 
planetary system. Most of the known planets in multiple planet sys- 
tems have masses roughly comparable to Jupiter (most likely due 
to detection biases). The unusual intermediate planetary masses re- 
vealed in the HD 69830 system have already inspired several the- 
oretical investigations. As a result, theorists have developed a de- 
tailed model for the mass growth and or bital migration of the ob- 
serve d planets through a series of papers ('Alibert et al.'200l l2006l : 
iLovi s et al. 2006). In particular, Alibert et al. ( 2006) found that the 
observed characteristics could be reproduced by starting with a gas 
disk in which the surface density, E, is related to the disk radius, 
a, by E oc tr^'- and normalized to 800g/cm^ at 5 AU (This sur- 
face density is around 4 times greater than the minimum mass solar 
nebula, giving Mrf,,,t = 0.07Mo (0.07 AU 30AU)). The disk has 
a dust-to-gas ratio of 1/70. They then inserted 3 planetary embryos 
of mass 0.6Me at initial semi-major axes of 3, 6.5 & 8 AU and used 
a semi-analytic model to migration and grwoth of the embryos due 
to their interaction with the disk over the course of the disk life- 
time, Tjisk - 2Myr. After accretion and migration, the 3 planets 
in the model were found to have shifted inwards to the observed 
semi-major axes, have total masses consistent with the minimum 
observed masses and possess core masses of ~ 10, ~ 7.5 & ~ lOMe 
respectively. 

While the I Alibert etai] ( I2OO6I) model does an impressive job 
of matching the observed planet masses and semi-major axes, it is 
a 1-D semi-analytic model, using numerous assumptions and ap- 
proximations to model the interaction between the protoplanetary 
disk and the growing and migrating planetary cores. This raises 
questions regarding the effects of the planets on any planetesimal 
disk and the feedback effects of disk evolution and "shepherding" 
on planetesimal accretion rates. In addition, the model does not ad- 
dress the eccentricity evolution of the planets and planetesimals. 
In this paper, we build on this work, as we continue the quest to 
understand the formation and orbital evolution of this fascinating 
system by applying n-body methods to follow planet esimal evolu- 
tion w ithin the context of the semi-analytical model of lAlibert et al.l 
( 12006) . 

1.2 IR Observations 

Spitzer observations of the HD 69830 system reveal a strong in- 
frared excess (relative to the stellar phot osphere) between 8 & 
35fim, but no significant excess at 70pra jSeichman et alj|2005h . 
This combination suggests that the system contains a disk of warm 
(~ 400K), small (< Ifiin) crystalline silicate grains orbiting close 
to the star (~ lAU). For these parameters, the coUisional timescale 
(for yum grains at ~ lAU from HD 69830) is ~ 4 00 yr, while the 
Poyn ting-Robertson drag timescale is ~ 700 yr l lBeichman et al.l 
l2005h . Thus, collisions in such a disk would grind down the grains 



until they became small enough to be removed via radiation pres- 
sure. Unless we happen to be observing the system at a very spe- 
cial time, this implies that the grains are replenished on a timescale 

< 1000 yrs so as to sus tain the observed IR exc ess. 

Several mo dels teeichman et al.l l2005l : luisse et alj |2007| : 
IWvatt et aill2007h have been proposed to explain the origin of the 
dust responsible for the IR excess, including: 

(i) A massive cometary population, 

(ii) The capture of a super-comet onto a circular orbit at ~ 1 AU, 

(iii) The steady-state evolution of a planetesimal belt at ~ 1 AU, 

(iv) A recent, large collision in a planetesimal belt at ~ 1 AU, 
and 

(v) Recent dynamical instability which results in planetesimals 
from an outer belt being thrown inwards. 

In the first scenario, the observed dust would be released by 
pristine comets enterin g the inner region of the planetary system. 
iBeichman et al. I ( l2005l) noted the observed spectral signature was 
similar to that of comet Hale-Bopp, but that reproducing the emis- 
sion using cometary ejecta similar to that of Hale-Bopp would re- 
quire ~ 10^ such comets per year to be delivered to the inner re- 
gions of the system. If the putative comets were of the same mass 
as Hale-Bopp and the dust were to be sustained for ~ 10' yr, then 
this would imply ~ 900Me of comets entering the inner regions of 
the planet system, unfeasibl y large for a residual Ku iper-Belt. 

In the second scenario, IBeichman et alj ( |2005|) suggested that 
a single object of the size of a large Kuiper Belt Object (e.g., Sedna) 
composed of ice and rock may have been scattered inwards to 
~ 0.5AU. They estimate that such an object would have an evapo- 
ration timescale of ~ 2 Myr, allowing a reasonable chance of ob- 
servation. Given the masses and orbits of the three known planets, 
a Sedna-like object at ~ 0.5AU could not be dynamically stable, 
unless it were on a low eccentricity orbit. In order to have a rea- 
sonable chance of observing the IR excess from such a body, the 
dynamics of the system must be such that a Sedna-like object could 
have recently been perturbed into the inner system (presumably on 
a highly eccentric orbit) an d then had it s orbit circularized. 

For the third scenario, IWvatt et alj j2007h considered the pos- 
sibility that the dust originates in a primordial planetesimal belt that 
is coincident with the dust and which has evolved in a quasi-steady 
state. They found that collisional processing would have removed 
most of the belt's mass over the > 2 Gyr age of the system and that 
the observed levels of dust are incompatible with this interpretation 
for a similar reason. One possible resolution to this is that plan- 
etesimals (and dust) may orbit with significant eccentricities, thus 
prolonging their survival, meaning that the emission at 1 AU could 
come from the inner edge of an extended disk when the planetesi- 
mals and dust are at pericentre (Wyatt et al, in preparation ). 

Similarly, for the fourth scenario [Wvatt et all ( |2007|) showed 
that an in situ planetesimal belt would be extremely unlikely to have 
undergone the single recent col lision that would gi ve rise to the ob- 
served IR excess. The results o f^ Lohneetal] ( l2008 l) may slightly in- 
crease the probability of a collision having occurred, as their model 
predicts higher remnant masses at late times, but a significant in- 
crease in probability would seem to require unfeasibly high initial 
disk masses. 

Finally, in the fifth scenario. IWvatt et al] ( |2007|) have proposed 
a model in which the current IR excess is due to a delayed epoch 
of orbital instability among planetesimals in an outer disk, perhaps 
similar to tha t of the Late Heavy Bombardment (LHB) in our own 
solar system lOomes et aT I ( I2OO5I) . Clearly this model then raises 
further issues such as how such a delayed instability might have 
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been triggered and why there has been no detection of the source 
population of planetesimals. 

All of these scenarios have unresolved problems, but their uni- 
fying feature is that they all depend to some extent on the existence 
of an extended distribution of planetesimals beyond the observed 
planets in the system. Fortunately the architecture of the planetary 
system sets constraints on how it formed, which in turn has impli- 
cations for the planetesimal population remaining following planet 
formation. By combining planetary formation models (described 
in gi.ll l. with the N-body techniques described in ^ we aim to 
develop an understanding of the remnant disk structure likely to 
populate the systems at the end of the planet formation process. 



1.3 Outline 

We perform dynamical simulations to investigate the formation and 
current state of the HD 69830 planetary system. Using numerical 
methods described in ^ we look in ^at the implications of the 
lAlibert et alj ( l2005h model for the planetary eccentricities and ac- 
cretion rates. In ^ we perform additional simulations to model 
the interaction between the planets and a putative protoplanetary 
distribution of planetesimals in order to model the predicted plan- 
etesimal distribution after the migration and growth of planets has 
been completed. In ^we consider possible variations in the forma- 
tion and orbital evolution of the HD 69830 planetary system. In ^ 
we consider the planetesimal distributions resulting from our sim- 
ulations and discuss their implications for explaining the observed 
dust emission. Finally, we summarize our conclusions and future 
prospects in 



2 METHODOLOGY 

We adopt the MERCURY n-body package of' ChambersI ( Il999h for 
the basis of our simulations, allowing us to follow the evolving or- 
bits of a population of planets and planetesimals. We model the 
planets as massive bodies, but treat the planetesimals as massless 
bodies (test-particles), adding in various additional routines to al- 
low the modelling of additional physical forces and effects such as 
i) Forced-Planet Migration, ii) Gas Drag on bodies of arbitrary size, 
iii) Dynamical Friction and iv) Mass Growth. All of these effects 
are implemented within MERCURY within the "mfo.user.for" sub- 
routine, thus ensuring that they are called sufficiently frequently 
that perturbations in mass and / or semi-major axes are adiabatic. 

To produce forced migration in a planetary system we imple- 
ment the req uired a i n ME RCURY using the constant-migration- 
rate model of IWvat such that the change in velocity, v is 
given by V = 0.5^ V GM* la^-j^ , where a: is a constant and M-,, is 
the stellar mass, giving a constant migration rate of 
a = 2\ ^[cFfGMl = k. This method allows us to arbitrarily select a 
constant migration rate or to adopt a time dependent migration rate 
using as an input the data from previous work, e.g. jAlibert et al.l 
l2006h . For the latter, we take a list of specified semi-major axes 
as a function of time and calculate from this list the resultant rate- 
of-change of semi-major axis, ii at any given time. This is then 
used as a time dependent k input to the afore mentioned migration 
model. We implement a similar mechanism to increase the plane- 
tary masses, where an M term can be specified and implemented as 
a function of time. 

Given that much of our work in volves migration rates dictated 
by the results o f lAlibertetai] ( l2006h . the planets are already subject 
to a force which gives rise to inward migration. We take this to 



specify the relevant timescale, and th en apply the simpl ified model 
for a, e and i damping described in IZhou et al.l ( l2007t) . in which 
any bodies of mass M orbiting in a gas disk background can be 
taken to suffer a gravitational tidal drag force, where to leading 
order in e & i, the semi-major axis evolution can then be linked to 
the eccentricity and inclination damping rates via: 

1 



Tjidal 

<e> 



5e2 + 2i-2 
1 



Tndal 

1 



We then implement the eccentricity and inclination damping in 
MERCURY in the mfojiser.for subroutine via direct damping of 
the relevant elements. 

We implement g as drag on the planet esimals using a model 
similar to that used bv lMandell et al.l ( l2007h . We consider a three- 
dimensional gas disk model where the gas density, p,,, and the ver- 
tical scale height, za, are given by (M andell et al..2007i') . 



0.0472 f—^ 
\\AU 



AU. 



(1) 

(2) 



where e = 11/4, y = 5/4 and the gas disk is taken to decay ex- 
ponentially on a tim escale Tjisk = 2Myr. Note that the model of 
lAlibert et alj ( |2006|) assumes that the 0.6M^ embryos will have 
taken ~ 0.93 Myr to grow. As such, our simulations are taken to 
start at f = 0.93 Myr, and hence any gas disk in the simulation will 
have dissipated by an amount e^" '^Z'"''"'. at the start of the n-body 
simulation. 

La rge bodies orbi ting in such a fluid will be subject to a decel- 
eration ( lRafikovll2004h 

d\ 3C PgVr 



dt An Pprp 

leading to a drag timescale of 



Tgd ~ 5 



PgVr 



(3) 



(4) 



where is the velocity of the particle w.r.t. the local gas velocity 
and the subscript, p, refers to the particle in question and we have 
assumed that planetesimals ^ 100km will have a drag coefficient 
C~\. 

The sub-keplerian local gas velocity, Vg, is calculated by con- 
sidering only the horizontal component of the stellar gravity in cal- 
culating the circular velocity and combining this with a further re- 
duction due to internal pressure support, giving 



(1 - '7)vk, 



where vk is the standard Keplerian circular velocity, 

, = (^)(|)(,_^^^)a„d/5 = 1/2 sets the dtsk temperature 

structure, T oc r"''. 

The above formalism makes it clear that a highly eccentric 
and/or inclined planetesimal will pass through regions of the disk 
significantly more removed from the centre and/or midplane of the 
gas disk, thus experiencing significantly lower gas densities and 
associated drag forces. 

Unless otherwise specified, we apply gas drag as though the 
planetesimals have a characteristic radius ~ \QQkm. We note that 
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this planetesimal size will cause the escape velocity from the 
planetesimals to be comparable to the velocity dispersion of the 
planetesimals for < 2i >x< e >~ 0.05, hence giving a quasi- 
equilibrium state jLufkin etaLll2006) . We take the initial distri- 
bution of our planetesimals to be a Rayleigh distribution with 
< 2i >=< e >= 0.05 

Whilst the addition of non-conservative forces to the MER- 
CURY routine prevents the traditional accuracy checks on conser- 
vation of energy and momentum, we did confirm that, (i) in the ab- 
sence of migration, using MERCURY in hybrid mode to simulate 
non-migratory planets embedded in a planetesimal swarm would 
typically result in fractional energy changes of ~ 10"'" and frac- 
tional ang. mom. changes of ~ lO"'"*, and (ii) when applying our 
gas damping routines to both non-migratory and migratory plan- 
ets in a static (non-dissipating) disk, we observed the damping 
timescales for a, e and i to scale as expected with mass, semi-major 
axis, eccentricity, inclination and gas density. 

Using the above algorithms in conjunction with the basic 
MERCURY package allows us to investigate a wide variety of 
physical effects in the HD69830 system. 

As a simplified test of the model, we conduct simulations sim- 
ilar to those of Lufldn^taL (2006), in which a Jupiter-mass planet 
migrates inwards at a rate of 10"^ AU per year through a planetesi- 
mal disk distributed between 0.5 & 4.5 AU. Using the rapid hybrid- 
symplectic (H-S) algorithm we typically find that the inward mi- 
gration of the planet will trap those planetesimals at smaller semi- 
major axes into mean motion resonances (MMRs) and then sweep 
them into more and more highly eccentric orbits, but that as the 
migration continues, many of the planetesimals are scattered out of 
the resonance into highly eccentric bound orbits, with an average 
eccentricity of 0.5. Repeating the simulation using the slower but 
more accurate Bulirsch-Stoer (B-S) algorithm results in a system 
in which the planetesimals remain tightly confined to resonance, 
being excited to higher and higher eccentricities as the planet shep- 
herds them inwards until finally, many of the planetesimals are 
ejected from the system, leaving fewer of the planetesimals surviv- 
ing in the system, and those that do survive have a lower average 
eccentricity (e ~ 0.1). 

The difference appears to stem from the way in which the 
two algorithms resolve the orbits of highly eccentric planetesimals 
trapped in the MMRs ahead of the planet. The B-S algorithm, re- 
solving the encounter in greater detail, results in much greater num- 
bers of planetesimals being confined within the resonance until they 
finally either (i) directly impact the planet, or (ii) are given a large 
kick and ejected from the system entirely. In contrast, the H-S al- 
gorithm tends to use a fixed timestep which is too large to properly 
resolve the fast moving, small pericentre planetesimals. This re- 
sults in many planetesimals "leaking" from the resonance, not be- 
ing shepherded to such high eccentricities, and subsequently tend- 
ing to only suffer high impact-parameter scattering and thus remain 
on bound (but excited) orbits within the system. 

W e note that the results originally reported by iLufkin et al.l 
( l2006l) are intermediate between the results of our B-S & symplec- 
tic tests, with few of the particles in their simulations being excited 
to the point of ejection. We thus suggest that the use of approximate 
integration techniques (such as the MERCURY H-S integrator and 
the PKDGRAV Runge-Kutta 4th order integrator) can, in cases of 
high eccentricity excitation, lead to anomalous results, and should 
be avoided when studying cases where significantly eccentric plan- 
etesimal populations arise. 

We also tested MERCU RY on simulations similar to those of 
lEdgar & Artvmowiczl j2004h in which a planet migrates from 6 to 



0.5 AU, scattering planetesimals initially distributed in an annulus 
between 3 and 4 AU. These simulations result in excitation of the 
planetesimals, but to a much lower level than in the previous test, 
with little or no ejection from the system occurring. We find that us- 
ing MERCURY in both H-S and pure B-S mode gives essentially 
identical results, and that these results follow exactly the same de- 
pendencies on planetary ma ss and migration rate as observed by 
lEdgar & A rtvmowiczl l l2004h (although we note that the absolute 
values of the eccentricity excitati ons we observe are always off -set 
slightly below the values found bv lEdgar & Artvmowiczl ( |2004|) us- 
ing the Runge-Kutte based PKDGRAV integrator). 

Having tested the integrators in this manner, we felt confident 
in the behaviour and accuracy of the H-S algorithm only when the 
planets are at large semi-major axes and have excited little eccen- 
tricity in the planetesimal population. For the purposes of speed, 
we therefore used the H-S algorithim when the planets are at large 
semi-major axes (when we were sure that the fixed timestep would 
be able to resolve the smallest pericentres) but then switch perma- 
nently to the B-S algorithm as soon as any of the planets migrate 
inside 1 AU and start to excite significant eccentricities and small 
pericentres. 



3 DYNAMICAL MODELLING OF PLANETARY 
ECCENTRICITY EXCITATION 

Here we consider a MERCURY model containing just 3 planets. 
3.1 Alibert Model 

We consider the initial conditions as adopted in the model of 
lAlibert et alj ( 12006) . We then force the three planetary embryos to 
grow in mass and decrease in semi-major axis according to the out- 
put of their model (kindly provided in numerical format by Yann 
Alibert). This allows us to assess the planetary eccentricities ex- 
cited by mutual perturbations during inward migration along the 
track proposed bv lAlibertet"ZI ( l2006l) . 

We conduct a large number of such simulations, each sim- 
ulation starting with the planets on orbits having zero eccentric- 
ity, small (< 1 degree) randomised inclinations and different ran- 
domised mean anomalies. A sample of some typical results is 
shown in Fig.[T] combining and averaging the results from 10 sim- 
ulations. As the planets migrate inwards, we find that the average 
planetary eccentricities excited along this forced migratory path are 
0.007, 0.008 and 0.013 respectively. We note that there is a large, 
fast increase in eccentricity visi ble in the m i ddle a nd and outer 
planets in Fig[T] In contrast with lAlibert et al.l j2006l) . we find that 
the crossing of the 3:2 MMR by the outer two planets is the source 
of this excitation. 

We note that the observed eccentricities and their l-cr errors 
were 0.10 + 0.04, 0.13±0.06 & 0.07±0.07 respectively, so the sim- 
ulation results are an order of magnitude below the mean observed 
values. Our own Markov-Chain Monte-Carlo (MCMC) analy sis of 
the radial velocity data using the methodology described in iFor j 
I2OO6) is summarised in Table[T] It again suggests that actual plan- 
etary eccentricities are likely to be significantly higher than those 
seen in our simulations. However, we caution that there can be sig- 
nificant uncertainties in the radial velocity observations. In particu- 
lar, best-fit orbital solutions can overestimate the eccentricity, e, of 
a nearly circular orbit, particularly for planets with small velocity 
amplitudes (Shen & Turner 2008). Zakamska et al. (in prep) find 
that the summary statistic e + P is significantly less biased 
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Table 1. MCMC Analysis of Planetary Eccentricities 



1 Planet 


1 




Percentiles 








2.5 


16.3 


50 


83.7 


97.5 


HD6983()-b 


0.010 


0.045 


0.088 


0.131 


0.173 


HD69830-C 


0.007 


0.042 


0.106 


0.174 


0.241 


HD6983()-d 


0.003 


0.020 


0.066 


0.141 


0.240 



than several other estimators, where h = e cos o), k = e sin ai, oj is 
the argument of periastron, and the tilde denotes the median value 
from a Bayesian posterior sample. Using the same MCMC simu- 
lations as reported in Table 1, we find = 0.076, = 0.080, 
and = 0.021, suggesting that the bias is modest for the inner 
two planets. In addition to the statistical errors reported in Table 1, 
eccentricities can be overestimated due to model misspecification 
(e.g., the presence of an additional planet that was not included in 
model fitting). Thus, we caution that further radial velocity obser- 
vations are warranted to define how non-circular the orbits are. 



3.2 More Massive System 

Given that the lAlibert et al.l ( l2006h model is calibrated to produce 
the minimum masses for the HD69830 system, any relative incli- 
nation between the system and the plane of the sky would mean 
that the actual planetary masses could be significantly above this 
minimum. Greater planetary masses during the migratory stage 
could allow the planets to self-excite to much greater eccentricities. 
Therefore, we run a new suite of simulations to look at the effect 
of scaling-up the initial embryo masses, as well as the subsequent 
planetary mass at each time-step in the simulations. For each mass 
scaling we run ten sets of simulations that differ only by the initial 
mean anomalies, giving an approximate measure of the variabil- 
ity of the potential evolutionary scenarios and outcomes. We stress 
that these increased mass mode ls are only an appr oximation (di- 
rectly increasing the mass in the lAlibert et alj j2006l) model would 
change various relative growth and migration timescales), but they 
do allow us to gain some insight into the behaviour of a more mas- 
sive system. 

We plot in Fig |2] the mean, upper-quartile and lower-quartile 
eccentricities from simulations for the three planets as we vary both 
the initial planetary masses and the subsequent mass growth rates. 
Over the top of this we also plot (solid lines) the observationally 
inferred eccentricities which result from our MCMC analysis of 
the best fit orbital elements to the radial velocity data. 

If we focus on the results for the inner two planets, then we 
can see that (i) values of e a; 0.01 (as seen in Fig[T) are approxi- 
mately coincident with the 2.5th percentile line, and (ii) the lower 
mass scalings (2x & 3x minimum mass) for both planets fall be- 
low the 16.3th percentile line. So, whilst it is not impossible that the 
observations are consistent with zero eccentricity, the likelihood is 
that it is significantly above this, suggesting that if one wishes to as- 
crib ethe observed planetary eccentricities to mutual planet-planet 
excitation during planetary migration, then the planets would be ex- 
pected to have masses over 5x the minimum mass from the Alibert 
models. Turning to the results for the third planet, and performing a 
similar analysis, we find that we expect from the MCMC fits a mass 
scaling less than 10 times the minimum mass. Taken together, these 
results suggests that we concentrate on mass scalings 5-10 times 
higher than the minimum mass models. Such a large increase in the 
masses would demand a nearly face-on system (4 < ; < 9 degrees 
for -j^ < sin i < ^), the probability of which is rather low, as well as 



the implication that the initial disk mass and surface density would 
be challengingly high. 

3.3 Gas Damping 

Given that the planets are growing and migrating within a gas disk, 
one should include gas-induced eccentricity damping in the system 
simulations. We model this using the planetary damping model de- 
scribed in ^and plot the results in Figure[3] Whilst the eccentrici- 
ties of the middle and outer planets are again clearly still excited by 
the 3:2 MMR resonance crossing at ~ 3 x 10^ yrs, this is damped 
back down on a timescale of a few x 10^ years. 

We find that this rate of damping is far too high to allow any 
significant eccentricity to be preserved in the system. For both the 
standard mass simulations (not shown) and for the simulations with 
5 times greater masses (Fig[3} we find that the eccentricities for all 
of the planets are damped to values lower than 10"^, two orders of 

magnitude below the observed values. 

Th ese results illustrate t he problem, e.g. IPapaloizou et al.l 

( l200ll) : IChatteriee et al.l J2007h . of how to generate an extra- solar 
planetary system in which there is sufficient gas present in the disk 
to cause mass growth and semi-major axis decay, whilst simulta- 
neously ensuring that the gas present in the system does not effec- 
tively circularise th e planetary orbits. Given the difficulties of the 
lAlibert et alj l |2006l) model, another form of eccentricity excitation 
must be posited following the decay of the gaseous disk. 



4 EXCITATION OF PLANETESIMAL ECCENTRICITIES 

There are a number of papers iFogg & NelsonI [20051. 12007 Jblf3 
lMandelletaIll2007h which focus on the formation of terrestrial 
planets during and after the migration of a Jupiter-mass planet 
through the inner system. They typically find that the passage of 
the giant planet does not clear the inner system of solids, but rather 
initially shepherds the material inwards before exciting it and fi- 
nally scattering/expelling it to exterior orbits. If this situation can 
be replicated in simulations of the HD69830 system, it could pro- 
vide a means of producing a scattered disk external to HD69830d 
which may result in a long-lived population of hot dust as observed. 

In this section we add massless planetesimals to out MER- 
CURY simulations of planetary migration. To investigate the scat- 
tering of planetesimals in this model, we perform a number of dif- 
ferent simulations to understand the effects of various physical phe- 
nomena on the final distribution of planetesimals after the conclu- 
sion of planetary migration. We vary the planetary masses, the gas 
drag acting on the planets and the gas drag on the planetesimals. 
The details of these simulation are recorded for convenience in Ta- 
ble H 

Unless otherwise stated, all simulations have an initial distri- 
bution of planetesimals which follows the minimum mass Solar 
nebular model surface density profile of E oc a"''^, with the ini- 
tial semi-major axes limited between 0.1 and 9 AU, and with the 
eccentricities and inclinations drawn from a Rayleigh distribution 
with 2 < i >=< e >= 0.05. 

In general, we perform composite simulations consisting of m 
different simulation runs, each simulation containing n planetesi- 
mals and 3 planets. This allows us to gain a speed benefit in finding 
the distribution of m x n planetesimals by running simulations in 
parallel, but more importantly allows us to include the range of dif- 
ferent planetary eccentricity excitations seen in 33.1| bv using ran- 
dom initial planetary mean anomalies in each of the parallel runs. 
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Figure 1. Eccentricity evolution for the three planets in the HD69830 system. Ten simulation runs are superimposed. The individual runs are plotted in grey, 
whilst the overall average is plotted in black. The planets grow and migrate according to the minimum mass Alibert model, with no gas damping operating. 
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Figure 2. The dashed lines and data points give the eccentricities for the thi'ee planets as functions of the mass-scaling (x-axis) in our simulations, each point 
gives the mean eccentricity over 10 identical runs, along with the upper and lower quartile bounds. The solid lines give the MCMC 2.5, 16.3, 50, 83.7 & 97.5 
percentile values for fits to the R.V. observations (50th percentile in bold, 16.3 & 83.7 in medium and 2.5 & 97.5th in light). 



thus giving us a view of the probabihties of different outcomes for 
a given set of starting conditions. 

4.1 Planetesimal Scattering in the 3-Planet Model 

Given the results of ^ we initially examine the case of planetes- 
imal excitation in an undamped system consisting of three planets 
following the growth and migration paths of 33.11 This is our sim- 
ulation set A. 

The shepherding effect on the planetesimals is clearly visible 
in the sample illustrations given in Fig [4] Here we see that as the 
planets migrate inwards, numerous planetesimals are forced to mi- 
grate ahead of the planets, subsequently rising up a number of the 
MMRs interior to the planetary semi-major axes, with the vast ma- 
jority being shepherded inside the innermost planet. From around 
0.5 Myr onwards, we start to see (Fig 4(g)) planetesimals being 
scattered out of the MMRs, some colliding with planets, some 
more being ejected from the system, whilst others are scattered onto 
high eccentricity orbits at intermediate semi-major axes. By ~ 0.8 
Myr, we find that the resonance structures interior to the innermost 
planet have now been almost completely destroyed, leaving a large 
swarm of planetesimals effectively populating the majority of the 
parameter-space 0.1 < a < 10, 0<e< 1. 

When the simulations have completed (I.e. the three planets 
are at the semi-major axes observed) we find that 27% of the plan- 
etesimals survive in the system, the rest being either ejected or suf- 
fering collisions with the planets or the central star. 



It should be noted that in the interests of brevity, these simu- 
lations were terminated after 10* years. If (impractically) the sim- 
ulations were continued on to billion-year time scales, then it is 
clear that many of the remaining planetesimals whose orbits cross 
those of the planets would also be ejected or suffer a collision, thus 
effectively clearing out most material with pericentre inside \AU, 
altho ugh some material may rem ain in a belt between p lanets c & d 
(Lo vis et alj2006l . |ji et alj2007h . However, as shown in lWvatt et afl 
cootI) . any material inside \AU which was sufficiently close to 
circularised to avoid planetary collision / excitation would be col- 
lisionally ground down by the current age of the system. As such, 
we therefore concentrate our subsequent analysis on those planetes- 
imals with pericentres outside 1 AU. 

We justify the neglect of the planetesimals with q < \AU hy 
now considering the possibility that the planetesimals discarded by 
this approach (i.e. the planetesimals which survive within the sys- 
tem at 1 Myr with q < \AU) could be scattered out to large semi- 
major axes / pericentres through interactions with the planets, thus 
enriching the extended eccentric disk. To investigate this scenario, 
we continue the integration of Set A from 1 Myr to 100 Myr, fo- 
cussing solely on the 12% of planetesimals which have ^ < 1 at 1 
Myr. We find that over the course of the subsequent 100 Myr in- 
tegration, over half of these planetesimals collide with the central 
star, a third collide with the various planets and a negligible per- 
centage are scattered out to the external disk. We summarise this 
data in the Table|2]under the entry labelled A+, where we have as- 
sumed that the amount of material with q > 1A{/ is fixed. Remov- 
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Figure 3. Eccentricity evolution for the three planets in the HD6983() system. Planets are 5x the minimum observed mass and are subject to gas damping. Ten 
simulation runs are superimposed, the individual runs plotted in grey, the overall average plotted in black. 



Table 2. Summary of the main planetesimal scattering simulations. 



Simulation Number of Planet Gas Damping % Ejected % Accreting onto % Surviving Surviving 

Planetesimals Mass on... or planets with with^>l 

(X Min.) Planets Planetesimals Hit Star b c d 1 < ^ % < e > < q> 



A 


10,000 


1 


No 


No 


4 


65 


3 


1 


12 


15 


0.30 


3.55 


A+ 










11 


66 


5 


2 


1 


15 






B 


2,500 


5 


No 


No 


43 


37 


1 


0.2 


3 


15 


0.47 


3.86 


CI 


2,500 


1 


Yes 


No 


4 


65 


4 


0.3 


12 


14 


0.29 


3.66 


Dl 


2,000 


1 


No 


Yes, 100km 





3 


1 


0.2 


80 


15 


0.13 


3.75 


D2 


2,000 


5 


No 


Yes, 100km 


12 


1 


0.5 


0.1 


74 


13 


0.28 


4.32 


El 


1,000 


1 


No 


Yes, 1000km 





6 


3 


0.6 


69 


23 


0.25 


2.92 


E2 


1,000 


5 


No 


Yes, 1000km 


34 


5 


0.6 


0.3 


41 


19 


0.46 


3.97 


F(i) 


2,300 


1 


No 


No 


17 


56 


2 


0.3 


11 


14 


0.29 


3.53 


F(ii) 


300 


1 


No 


No 


18 


57 


2 





9 


14 


0.29 


3.47 


F(iii) 


700 


1 


No 


No 


24 


49 


3 





11 


13 


0.27 


3.97 


F(iv) 


100 


1 


No 


No 


44 


34 


1 


1 


9 


11 


0.30 


4.20 


G 


2,500 


1 


No 


No 


40 


38 


2 





9 


11 


0.41 


3.83 


H 


1,250 


1 


No 


Yes, 100km 





0.3 


0.8 


1.2 


84 


13 


0.15 


3.59 



ing those planetesimals whose pericentres are inside 1.0 AU, leaves 
us with 15% of the original planetesimal material occupying orbits 
with pericentres greater than 1.0 AU. Note that the initial propor- 
tion of bodies in this same region was approximately 21%, so this 
region is depleted only slightly from its initial value. However, the 
mean eccentricity has been greatly excited to < e >= 0.30. See Fig 
|5(a)| for the eccentricity and pericentre histograms. We note that the 
number of planetesimals scattered to semi-major axes larger than 
10 AU is very small. 

We should note that the migration rate in the Alibert model 
is artificially redu ced compared to the analytic estimates of 
iTanaka et alj bOOA . a reduction of 1-2 orders of magnitude be- 
ing necessary to explain observations of extr a solar giant planet s 
(E.g.'Alibert et al.' ('2005');'Daisaka et al.' (l 2006h : llda & Lii] j2008h : 
[Senz et al. (2008)), and also suggested by various other results in 
which the migration behaviour is found to be more complex than 
that suggested by linear typ e-I theory. In particular, the results of 
iNelson & Papaloizoul ( |2004) suggest that migration can become 
stochastic as a result of disk turbulence. If such an instantaneously 



fast, stochastic migration was overlaid on a much longer timescale 
orbital decay, then the overall migration rate could be similar to 
that used in the Alibert model. However, the instantaneously fast 
stochastic migration would work to reduce shepherding, since the 
planetesimal trapping probability is a strong function of mi gration 
rate (w ith lower probabilities for faster migration rates, e.g.. lWvatll 
( l2003l) ). and stochastic variations could c ause previously trapped 
planetesimals to fall out of resonance (e.g.. Murrav-Cl av & Chiaii3 
(,2006.) : . Adams et al. c200&) ; ,Rein & Papaloizou (200^1 1 



4.1.1 Surviving Planetesimals 

The planetesimals which "survive" in the system are those plan- 
etesimals which are neither ejected from the system, nor suffer col- 
lisions with either the central star or one of the migrating planets 
(The planetesimals are treated as test particles so no planetesimal- 
planetesimal collisions take place). To gain some additional insight 
into the behaviour of the system, we look in detail at the planetesi- 
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Figure 4. Evolution of eccentricity - semi-major axis for tlie planetesimals distribution due to tlie forced migration of planets. Results of Simulation Set A ■ 
Minimum mass planets with no gas damping present in the system. 
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(a) Simulation A: Standard Alibert-Mass Simulations after 1 Myr 
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(d) Simulations D & E - All have Standard Mass Planets. LEFT - Eccentricity v. Semi-major Axis for planets and planetesimals for Set D only. CENTRE - 
Eccentricity Histograms for Gas Drag on 100km Planetesimals (Set D, Solid Histogram), Gas Drag on 1000km Planetesimals (Set E, Dashed Histogram) and 
Standard Zero Drag (Set A, Dotted Line). RIGHT - Perieentre Histograms, Line-types as for centre plot. 

Figure 5. Comparative results plots for selected simulations. Left Hand Column = Eccentricity-v-Semi-Major Axis Plots for Planets & Planetesimals; Central 
Column = Final Eccentricity Histograms for planetesimals (only showing planetesimals with q < 1); Right Hand Column = Final Perieentre Histograms for 
planetesimals (only showing planetesimals with q < I). 



mals that do not survive to find out if and how their ultimate fate is 
dictated by their initial starting positions. 

We plot the fate of the planetesimals as a function of their 
starting semi-major axis in Fig|6] This makes it clear that the plan- 
etesimals which are located well within the initial orbit of the inner 
planet (3 AU) have a higher chance of being lost from the system 



(~ 70% for the inner-most planetesimals), whilst those that are lo- 
cated outside 2AU are almost certain to survive as bound bodies 
within the system. In addition, we see that almost all of the plan- 
etesimals which start outside 2 AU will end up on an eccentric orbit 
with a pericenter beyond ~ lAU. Thus, these planetesimals could 
be "useful" for providing a parent population that explains the ori- 



10 M.J. Payne 




0.1 1 10 



Start Position / AU 

Figure 6. Fate of planetesimals as functions of starting semi-major axes 
for: (i) Fraction of planetesimals which survive with q > 1.0 - Solid Line; 
(ii) Fraction of all planetesimals which do not survive - Dashed Line; (iii) 
Fraction of planetesimals which do not survive becasue they hit the inner 
planet - Dotted Line. Excludes any gas damping effects. 

gin of the 3 - 35yum emission (contingent on their eccentricities and 
collisional lifetimes ). 

Next, we investigate in detail the fate of the individual plan- 
etesimals. We find that out of the 10,000 planetesimals with which 
we start the simulation, 65% hit the inner planet, 3% hit the middle 
planet and 0.5% hit the outer planet. A further 3% were put into 
parabolic orbits such that they collided with the central star or were 
ejected from the system. We note from Fig 4(i) that the majority 
of the collisions with the inner planet are happening between 0.7 
& 0.8 Myr into the simulation, at which point the inner planet is 
effectively pushing into and colliding with the planetesimals in the 
MMRs ahead of it. To eliminate the possibility that these results 
are unduly influenced by the surface density profile that we have 
adopted in the simulations, we repeat the analysis for a set of sim- 
ulations in which the surface density profile is distributed evenly 
between 0.5 and 9.5 AU and in which we try removing 1 or 2 of 
the outer planets. These simulations showed that the surface den- 
sity profile does not skew the results. We find a very similar ejection 
profile to that with a E oc a"^'^ profile, with the planetesimals which 
do not survive as free bodies within the system primarily suffering 
collisions with the inner planet. 

4.1.2 Accretion Rate 

Given that the solid disk in the Alibert model is massive, containing 
~ 0.5Mj of solid material in the 0.1-9 AU region that we simulate, 
an impact fraction of 65% onto the inner planet would imply that 
the solid core of the planet would attain a mass of ~ 0.3Mj or 
~ lOOMg,, far above both the observed minimum masses for the 
system and the ~ lOMe cores in the Alibert model. We compare 
the growth rates for the planets as inf erred from the plane tesimal 
accretion rates with those found in the lAliberte7^ ( l2006h model 
in Figl?] We find that the inner planet, starting out at 3 AU would 
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Figure 7. Comparison of planetary growth rates inferred from standard ac- 
cretion scenari o (Set A, thin , soUd lines) with those from the semi-analytic 
growth rates of l Alibert et aljj2006l) (Thick, dotted lines). Black plots - Inner 
Planet(b), Red plots - Middle Planet(c), Blue plots - Outer Planet(d). 

collide with many more planetesimals than the outer two planets, 
thus growing 10 - 100 times more massive, significantly at odds 
with the observed mass ratios . 

The lAlibert et al.l ( l2005h model is a semi-analytic model in 
which the core accretion is taken to occur smoothly from an annu- 
lar region extending out to 4 Hill radii either side of the planet, with 
the disk profile depleting self-consistently, but always maintaining 
a local surface density profile oc tr'/-. In contrast, our dynamical 
model makes no such assumptions, simply recording the number 
of collisions between the planetesimal population and the planets. 
Fundamentally, we would expect this to give a more accurate rep- 
resentation of the accretion rates, as long as no crucial physics have 
been omitted. 

To reconcile the results of our n-body simulations with the 
semi-analytic results would require that the relative number of 
planetesimal collisions onto each of the planets in the dynamical 
model be in approximately the same ratios as the core masses of 
the three planets in the growth model. Given the ubiquitous na- 
ture of the shepherding phenomena in n-body simulations, it seems 
highl y unlikely that a ny simple semi-analytic model (such as that 
of .Alibert et al. Illool))' which neglects this fundamental rearrange- 
ment of the solid disk profile can ever be in agreement with n-body 
simulations, suggesting the necessity of making improvements to 
these semi-analytic models, in particular the treatment of shep- 
herded low eccentricity planetesimals. 

Finally, we note that the accretion rates in our model were 
calculated directly from the collision rate between the planets and 
planetesimals, and that this required that we make an assump- 
tion about the planetary density {r oc M''^p"'''). We used a con- 
stant density of i.5g cmr~ throughout our entire calculation, corre- 
sponding approximately to that observed in the Solar-System ice- 
giants, initially suggesting that our accretion rates at the start of 
the simulation (when the body is a solid core) will be slightly too 
high, while at the end of t he simulations they should be approx- 
imately correct. However, iFortier et al.l ( l2007h show that the ef- 
fect of gas drag in the planetary envelope increases mass growth 
rates by up to a factor of 2 compared to their standard assump- 
tion of a pure solid core with density 3.2g cnr~. I.e. they require 
a higher effective radius, or alternatively a lower effective density, 
peff = 3.2 X (j^ gcm"^ X l.lg cm"'. Thus the initial growth 
rate in our simulations may be more accurate than initially implied. 
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However, we emphasise that these shght uncertancies in the accre- 
tion rates would apply to all planets in our simulations and would 
not therefore help to reconcile the order of magnitude difference in 
accretion rates that we observe between the inner and outer planets. 



4.2 Planetesimal Eccentricity Excitation as a function of 
Planetary Mass & Planetary Eccentricity 



4.2.1 Planetary Mass 



We demonstrated in 93.2l that the mass of the planets can have a sig- 
nificant impact on the eccentricity that the planets excite amongst 
themselves as they migrate inwards. We now look at how the mass 
of the planets affects the properties of the scattered planetesimal 
disk. In simulation B, we look at a model in which the initial masses 
of the planets and their subsequent growth rates are 5 times those 
of the standard Alibert model. 

In Fig |5(b)| we can compare the distributions of bound plan- 
etesimals with pericentres greater than 1.0 AU for simulations A 
& B (where the masses are respectively 1 and 5x those of the stan- 
dard model). We find that as the mass of the planets is increased, the 
available fraction of planetesimals remaining bound remains static 
at 15%, but the eccentricity distribution is shifted, with the higher 
mass planets creating more higher eccentricity planetesimals at the 
expense of low eccentricity planetesimals, increasing the average 
eccentricity of this population from 0.30 to 0.47. The pericentre dis- 
tribution is not as greatly affected, although there is again a slight 
increase in the number of high q planetesimals at the expense of the 
number of low q planetesimals. In addition we find that the range 
of initial semi-major axes which is ejected is not changed by an in- 
crease in the planetary masses: again, only the planetesimals inside 
2 AU are ejected/suffer collisions with any efficiency, whereas out- 
side this region, almost all the planetesimals remain bound to the 
system as free bodies. 

4.2.2 Planetary Eccentricity 



We saw in 933] that gas damping significantly reduces the plane- 
tary eccentricities. When we look at the effect on the planetesimals 
though, the difference made by damping the planetary eccentrici- 
ties is nowhere near as significant: we see in Fig |5(c)| the results of 
simulation C 1 in which the planets are damped (but the planetesi- 
mals remain undamped). We find that the the distribution of plan- 
etesimals is very similar to that found from simulation A. Overall 
this suggests that the eccentricity of the planets has only a small 
influence on the final distribution of the planetesimals, probably 
because the vast majority of the excitation is caused by the inward- 
shepherding and excitation, and this takes place undisturbed, irre- 
spective of (modest) planetary eccentricity. 

4.3 Gas Damping of Planetesimal Eccentricity 

We simulate gas damping on the planetesimal population using the 
model described in ^ We again note that the start point of these 
simulations with the embryo mass at 0.6Me is taken to correspond 
to a "system-time" of 0.93Myr, i.e. the gas disk has already suffered 
a significant amount of dissipation. In addition we note that the dis- 
sipation of the disk and the growth of the planets is not calculated 
self-consistently. 

The Alibert model assumes a disk lifetime of 2 x 10* years, 
so we adopt this as our disk dissipation timescale. Since we adopt 
an exponential model, dissipation will result in a significant low 
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Figure 8. Eccentricity vs. Semi-Major Axis plot after planetary migration 
for damped planetesimals: BLACK - 100km (Set Dl); RED - 1,000km (Set 
El). 



density "tail" of gas remaining in the system and thus cause non- 
trivial damping for a few times this timescale. As such we run our 
simulations for 2x 10' years to allow the planetesimals to approach 
an equilibrium. 



4. 3. 1 Scattered Disk 

In simulation Dl, we apply drag as if the planetesimals had a phys- 
ical radius of 100km and a density of Ig cm"-, and assume that the 
planets have the standard minimum masses and are unaffected by 
gas drag. We find in Fig |5(d)| that we still have around 15% of the 
planetesimals remaining in the useful zone out beyond 1 AU but 
that this population now has the much lower average eccentricity 
of 0.13. At this stage in the simulation of the system's evolution, 
the gas density has decreased to ~ 10"^ of its initial value, hence 
the possibility of any further eccentricity damping can essentially 
be neglected. In simulation D2 (not plotted) we repeat this gas drag 
simulation for planets with 5x the standard Alibert masses, and 
now find that even with this significant gas drag present, the average 
eccentricities of the remaining planetesimals can have < e >~ 0.28. 

In simulations El & E2, we repeat the gas drag simulations 
of Dl & D2, but we now model the gas drag on the planetesimals 
as appropriate for bodies with a physical radius of 1,000km. From 
Eqn|4]we see that the gas drag timescale, 

i"GD 0^ rp, 

so the damping timescale for a 1,000km planetesimal will be over 
twice as long as for a 100km one. The overplotted results in Fig 
1 5(d) I and Fig [8] show that these larger bodies, less affected by gas 
drag, have a mean eccentricity of < e >= 0.25 after the migra- 
tion of standard mass planets (El). More massive planets (E2, not 
plotted) have < e >= 0.46. So the much larger planetesimals of 
radius ~ l,OOOA:m will tend to occupy orbits of almost double the 
eccentricity of the smaller bodies of radius ~ lOOkm. 

We emphasise again the generic implication of the above: Any 
initially mixed population of planetesimals subject to gas-drag will 
become size-sorted over time, with the larger bodies tending to oc- 
cupy the high eccentricity orbits, whilst the smaller bodies (Dl & 
D2) will tend to populate more circular orbits. 
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Figure 9. Comparison of planetary growth rates in gas-drag simulations 
(Set Dl, thin , solid lines) with those from the semi-analytic growth rates of 
lAlibert et all )20Q6h (Thick, dotted lines). Black plots - Inner Planet(b), Red 
plots - Middle Planet(c), Blue plots - Outer Planet(d). 



Table 3. Results of Simulations with an Additional Embryo Present 



Outcome 


L CI UCllLagC 

of 

Simulations 


Planetary 
Eccentricity 


Planetesimal 
Simulations? 


Collision with (b) 


10 


0.17 


No 


Collision with (c) 


61 


0.11 


F(i) 


Collision with (d) 


10 


0.27 


F(ii) 


Scatter to Outer System 


4 


0.23 


F(iii) 


Ejection from System 


9 


0.19 


F(iv) 


Other (Multiple Collisions) 


6 


0.45 


No 



In summary, we find that the eccentricity damping of planetes- 
imals does not act as an insurmountable obstacle to the existence of 
an excited eccentric disk with pericentres ^ 1 AU in the HD69830 
system. The amount of available material with pericentre > 1 AU 
remains essentially unchanged at 15%, but the average eccentricity 
is reduced to 0.13. This still means that we have more than 10% of 
the solid disk mass available in an excited form to potentially act as 
the source of the emission. Further work will be required to un- 
derstand whether the observed distribution can survive collisional 
processing for the > 2 Gyr age of the system (see ^for further 
discussion). 



and mo derately eccentric pl anetesimal disks. The first two indicate 
that the lAlibert et alj dlOOSl) model will need tweaking to provide 
a self-consistent formation model. The third provides significant 
grounds for optimism in producing a parent population of eccen- 
tric planetesimals, but as we do not yet know whether the degree 
of eccentricity found is sufficient to extend the collisional lifetime 
of the system past 2Gyr, we wish to take this opportunity to under- 
stand whether further eccentricity excitation can be fostered. Be- 
low we consider three ways in which the model may be tweaked 
and consider its implications for planetary eccentricities, planetary 
accretion rates and planetesimal eccentricity excitation. 



4.3.2 Accretion Rates 

In 34. l.ll we saw that the collision rate of planetesimals onto the in- 
ner plane t was far too high to be compatible with the semi-analytic 
model of ' Alibert et al] j2006l) . In E l.ll we saw collision rates for 
the three planets of 65, 3 & 0.5 % respectively, whereas with gas 
damping acting on the planetesimals the problem was somewhat 
ameliorated: collision rates reduced to 3, 1 and 0.2 % respectively. 

The difference is essentially due to the fact that in the damped 
case, all of the planetesimals in the MMRs inside of the inner planet 
are damped to e a 0, meaning that their orbits do not overlap with 
that of the inner planet and hence collisions rarely occur. However, 
even with gas damping on planetesimals, the number of planetes- 
imal collisions onto the inner planet is still an order of magnitude 
greater than onto the outer planet. 

The growth profile plotted in Fig|9]makes it clear that we still 
have an order-of-magnitude difi"erence in planetary masses (accre- 
tion rates), but now all of the planets are too low in mass (core 
masses ~ lOM^ would require accretion of ~ 6% of the solid disk). 



5 OUTSTANDING PROBLEMS AND POTENTIAL 
SOLUTIONS 

The work in sections [3] and |4] has highlighted a number of results 
and outstanding problems. These include, (i) Lack of excitation 
and possible over-damping of the planetary eccentricities, (ii) The 
incompatibility of the collision rates with the core masses from 
lAlibert et alj ( I2OO5I) and (iii) The existence of relatively massive 



5.1 Additional Embryos 

The large initial separation of the three embryos considered thus 
far and their relatively low masses would suggest that a nascent 
system in this configuration would have a large number of simi- 
lar size embryos distributed between them. Whilst it is possible to 
conceive of this intervening material being accreted in an adiabatic 
manner, it also seems reasonable to suggest that these "additional" 
bodies could be growing and migrating in parallel with the three 
bodies considered thus far, opening the door to catastrophic impact 
and / or ejection scenarios. We briefly investigate this scenario by 
conducting 100 simulations in which an additional planetary em- 
bryo of mass O.6M0 is initially randomly placed between 4. 1 and 

5.2 AU (i.e. in th e middle third of the gap between HD69830-b and 
-c). The standard [Alibert et alj ( 1200^ semi-major axis evolution is 
applied to the 3 standard planets, whilst a scaled Type-I migration 
rate (lOx slower than analytic calculations suggest) is applied to 
the fourth planet. This gives the additional planet a migration rate 
intermediate between that of HD69830-b & -c. Gas damping is not 
applied to the planets. 

We typically find (in 81% of simulations - See Table|3j that the 
fourth planet collides with one of the other planets in the system. 
This tends to result in rather excited systems with mean eccentric- 
ity > 0. 1, rather than the under-excitement that we have seen previ- 
ously. Reducing the mass of the added planet would, of course, re- 
duce the level of additional excitation given to the remaining three 
planets. 

For a limited number of the simulations, we perform further 
detailed simulations to calculate the expected planetesimal distri- 
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bution resulting from the addition of an additional embryo. These 
additional simulations, referred to as F(i)-F(iv), are detailed in Ta- 
bleland a summary of the results is plotted in Figure [Tol 

Irrespective of whether the additional planet collided with one 
of the other planets (simulations F(i) & F(ii)), is scattered to the 
outer system (F(iii)) or is ejected from the system (F(iv)), the pop- 
ulation of scattered planetesimals with g > 1 is very similar to that 
observed in the standard simulation (set A), having very similar 
population fractions and very similar eccentricity excitations. N.B. 
There are fewer particles in Sim F(iii) than in F(i), hence the differ- 
ence in appearances of plots 1 1 0(a)l and [TO(b) [ despite both having 
similar fractions of retained particles with q > I. The only signifi- 
cant difference appears to be the slight reduction in the number of 
planetesimals accreting onto the inner planet, the difference essen- 
tially being made up by an increase in the number of planetesimals 
ejected onto parabolic orbits. However, the accretion rate onto the 
inner planet is still an order of magnitude larger than that onto the 
other two planets. 

If we repeat the insertion of an additional planet, but now ap- 
ply gas damping to the planets, then the situation becomes very 
different. We now find that in 98% of cases the additional planet 
remains securely trapped between the two planets it starts between, 
migrating in between them, with any mutual eccentricity excita- 
tion being quickly damped away (average planetary eccentricities 
drop to < e >~ 0.0005). Whilst this scenario clearly will do little 
or nothing directly to stir or excite the outer population of plan- 
etesimals, it may provide a simple way to initially trap planetary 
embryos in the inner system before subsequent perturbations eject 
them into the outer system, potentially providing a trigger for an 
LHB-type scenario. However, without long term simulations this 
remains speculative, especially given that the low planetary eccen- 
tricities (~ 10"^) may preclude any significant secular evolution to 
drive an instability. 



5.2 Migration after Growtli 

If the planets carry out the majority of their mass growth whilst they 
are stationary, then the eccentricity that they self-excite as they mi- 
grate may be significant, perhaps even more than that found in sec- 
tion [3]2l as a significant fraction of the eccentricity growth occurs 
as a result of MMR crossing, and as such is dependent on plane- 
tary mass at the time of resonance crossing. To investigate this in 
simulation set G, we take the extreme case of the planets growing 
to full mass at their initial locations and then migrating through the 
planetesimal disk to their final location. 

We find in Fig [TO] and Table [2] that the distribution of the re- 
maining planetesimals is such that the percentage of planetesimals 
with pericentres > 1 AU is now somewhat reduced compared to 
set A (11% versus 15%), but the average eccentricity is now sig- 
nificantly higher (0.4 versus 0.25). Thus, as might be expected, the 
results are somewhat similar to those simulations with increased 
mass and growth rates (E.g. set B, D2 & E2) in that the average 
eccentricity has been increased, but this time we also have more 
efficient clearing of the disk external to 1 AU. 

Finally, we note that the high initial masses in this model serve 
to excite significant eccentricities as the planets migrate, resulting 
in average eccentricities (over 100 simulations) of 0.24, 0.17 & 0.09 
for the inner, middle and outer planets respectively. 



Table 4. Details of simulation outcomes for various different initial plane- 
tary semi-major axes 



Initial Semi-Major Axis 


Accretion Rate onto Planet 


Average 


(AU) 








(%) 


Planetary 


b c 


d 


b 


c 


d 


Eccentricity 


3.0 6.3 


8.0 


2.7 


1.0 


0.2 


0.001 


3.0 3.5 


4.0 


1.2 


0.2 


0.7 


0.22 


5.0 6.0 


7.0 


1.2 


0.2 


0.2 


0.23 


1.0 6.3 


8.0 


0.6 


1.6 


0.2 


0.09 


0.5 6.3 


8.0 


3.0 


2.3 


0.4 


0.09 


1.0 2.0 


8.0 


0.5 


1.1 


1.3 


0.0002 


0.5 2.0 


8.0 


0.9 


1.1 


0.3 


0.0001 


1.0 4.0 


6.0 


0.2 


2.0 


0.3 


0.05 


0.5 1.0 


6.0 


0.6 


0.8 


1.0 


0.0002 



5.3 Different Initial Semi-Major Axes 

Given that the inner planet suffers a much higher collision rate than 
the other two planets, we consider a scenario in which the inner 
planet starts closer to the central star. This means that it migrates 
through a much smaller portion of the disk and thus has the poten- 
tial to interact and collide with a smaller number of planetesimals. 
In addition, a greater number of planetesimals remain in the un- 
perturbed region exterior to the inner planet, thus increasing the 
chance of impacts occurring with the two outer planets. We keep 
the migration timescale the same as well as M{t), but alter a{t) such 
that the planets migrate inwards at a constant rate, a(t) = k, arriving 
at their final semi-major axes after 0.9Myr. 

When the effects of gas drag on 1 00km planetesimals are ex- 
cluded, rearranging the initial semi major axes of the planets does 
not have any significant effect, unless the inner planet is started 
from exceedingly small semi-major axes (due to the power-law na- 
ture of the planetesimal disk). When gas drag is included, then con- 
ducting a variety of simulations in which the initial positions of the 
three planets are varied gives the results shown in Table |4] The 
first three columns give the initial semi-major axes, the next three 
columns give the resultant accretion rates onto the planets whilst 
the final column gives the overall average planetary eccentricity. 
For reference, the top line includes the results from 33.31 where gas 
drag was included on the standard distribution of planets. The most 
promising results (with regards to the accretion rates) from the ad- 
ditional simulations seem to occur when the inner two planets (b 
& c) are started close together and a relatively small semi-major 
axes, whilst the outer planet starts much further out: see for exam- 
ple the last line of Table|4l in which the planets start at 0.5, 1.0 & 
6.0 AU respectively and the subsequent fractions of planetesimals 
accreting onto the planets are 0.6, 0.8 & 1.0 respectively. 

We note in passing that the eccentricities of the planets in these 
models (right hand column of TableO are rather harder to tie in to 
the observed values. The final row of the table (which has the most 
even distribution of planetesimal collisions) has planetary eccen- 
tricities which are around four orders of magnitude lower than the 
mean observed values. Clearly one of the other models, in which 
the outer two planets start closer together and hence pass through 
a significant MMR and thus generate eccentricities of ~ 0. 1 would 
be favoured if future RV measurements do confirm the current ec- 
centricities. 



14 M. J. Payne 




10 10,0 0.2 0.4 0.6 0.8 1 2 4 6 8 10 

Log (Semi-mojor o.is / AU) Eccenlrioity Pericenire 



(a) Simulation F(i): Additional Planet, Collides with HD69830(c) 




(c) Simulation G: Migration After Growth 



0.98 Myrs 









A. dotted: f=0.15, <e>=0.30 
Dl. solid: f=0.15. <e>=0.13 
H, dashed: r=0.13. <e>=0.15 



A. dotted: f=0.15, <q>= 3.55 
Dl. solid: f=0.15. <q>= 3.75 
H, dashed: r=0.11. <q>= 3.59 




Log {Semi-major axis / AU) 



(d) Simulation H: Comparison of simulation H (Gas Drag in operation + Initial Semi-Major Axes altered to 0.5, 1.0 & 6.0 AU) with the standard simulation, 

l^gure^lf. *S)i^^ai%tive resuFts ^ots for selected simulations. Left Hand Column = Eccentricity- v-Semi-Major Axis Plots for Planets & Planetesimals; Central 
Column = Final Eccentricity Histograms for planetesimals (only showing planetesimals with ^ < 1); Right Hand Column = Final Pericentre Histograms for 
planetesimals (only showing planetesimals with ^ < 1). 



6 DISCUSSION 

6.1 Planet Formation Models 

Numerous studies jKokubo et al.l2006l : lRavmond et alj200i) show 
that the growth of embryos to isolation mass in a S oc a"''^ disk 



will occur inside-out, that is, the timescale for growth to isolation 
in such a disk will be t, ;o oc and hence isolation will occur quick- 
est in the inner disk. Therefore a migrating embryo which initially 
formed at 3 AU will be migrating through a region of space which 
is not occupied by a smooth solid disk, but rather by a series of 
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massive embryos which have already grown to isolation. The core 
growth of the embryo must then take place via giant impacts, which 
will clearly have a significant impact on both the planetary semi- 
major axes and eccentricities which would be predicted at the end 
of the growth and migration regime. 

In addition, we have seen in our work above that significant 
"shepherding" of material can occur in front of migrating proto- 
planets. This phenomena combined with the stirring up of the plan- 
etesimal disk by the passage of the first embryo means that the core 
growth rate for the first embryo is (at least) an order of magnitude 
greater than for the subsequent cores, a result that is significantly at 
odds with the assumptions of semi-analytic models. 

We also find that embryos which have grown more massive 
prior to passage through MMR are prefered in the model, as this 
allows an easier, more natural explanation of the observed plane- 
tary eccentricities, although the subsequent gas damping of these 
eccentricities still presents problems. 

These results serve to emphasise that a consistent model for 
planet formation and evolution is not possible using the current 
generation of semi-analytic models. In order to make progress, the 
semi-analytic growth models need to be combined with n-body 
models to ensure that both the mutual planet-planet interactions 
as well as the planet-planetesimal interactions are self-consistently 
modelled and able to be compared with observations. 



tion of the planetesimals is extremely important, since it is only 
> IQOkm planetesimals that remain on highly eccentric orbits fol- 
lowing scattering by planet migration and subsequent gas damp- 
ing. Thus, it is important to consider the fraction of the remain- 
ing 25 - 45Me of planetesimals that have radii > lOOkm. The 
planetesimals would be expected to have a bimodal populat ion of 
the type seen in the simulations of iMorishima et alj ( l2008h . who 
found that at lO^yrs, ~ 70% of the solid mass would be in the 
form of a few large embryos, with the rest residing in a popula- 
tion of planetesimals with a number density distribution given by 
n(m)dm = rrT^dm. Such a distribution has mass evenly spread at 
all sizes (on a log-scale), and thus we would anticipate that ap- 
proximately 10% (~ 2.5 - 4.5Me) of the surviving planetesmal 
material would reside in bodies of size 100 - 1000 km, with a fur- 
ther 70%_(2_J^-32Mj5) locked up in the very large embryos. In 
the I Alibert et all ( 1200^ model, from the initial Q.SMj = 160Me 
of solid material in the disk, only ~ 27 Me (15%) is accreted onto 
the migrating cores, so the contribution of the large embryos to the 
scattered distribution could be significant. Thus there is likely to be 
a significant mass fraction in large scattered bodies, which, along 
with the size sorting of the planetesimals via gas drag, works to en- 
sure that a large proportion of the remnant planetesimal mass would 
be locked up in large bodies which would be on the most eccentric 
orbits. 



6.2 Eccentric Planetesimal Disk 

We find in ^that the end product of our n-body simulations is 
the prediction that the excitation and retention of an eccentric plan- 
etesimal disk exterior to the three planets in HD 69830 is almost 
certain to take place (within the model parameters investigated). 
Significantly, we found that the mass of planetesimals remaining in 
this disk was large: after the planets had migrated and the gas disk 
dissipated, there would be ~ 25 - 45Mgj of material present in the 
scattered disk (using the initial disk profile of ( Alibert et al. 2006)). 

The presence of this amount of material in the young system 
(few Myr) does not necessarily guarantee that significant amounts 
of mass will remain in the system at late times (few Gyr), as the 
steady-state collisional processing of the disk will remove a large 
amount of material. The results of IWvatt et al.l ilQQ l) indicate that 
low eccentricity planetesimal belts evolve to a mass that is inde- 
pendent of initial mass, and in the case of HD69830, if this mate- 
rial were placed in a ring at 1 AU, then after 2Gyr the collisional 
processing would result in just 10"' Me of planetesimals remaining, 
which would be 1, 000 times less than required to explain the ob- 
served emission. This conclus i on ma y be mitigated to some extent 
by the results of iLohne et alj ( 1200 Sl) which show a small depen- 
dence of dust mass at late times on initial planetesimal mass, but 
this would probably not be sufficient as there results imply that the 
initial planetesimal mass would have to be unrealistically large to 
explain the observed emission in HD 69830. 

To work out how our planetesimals evolve we need to con- 
sider the collisional evolution of a population with a range of ec- 
centricities and semimajor axes: This has yet to be achieved, how- 
ever, work in preparation (Wyatt et al) already shows that increas- 
ing the eccentricity of a planetesimal belt from nearly circular to 
0.9 while keeping the pericentre constant increases the amount of 
material remaining at late times by 1 to 2 orders of magnitude. Al- 
though certainly giving no guarantee, this increases confidence that 
the distributions we derived might be able to explain the observed 
emission. 

However, our simulations also show that the size distribu- 



7 SUMMARY 

We used the model o f lAlibert et all ( |2006|) as the basis for under- 
taking a suite of n-body simulations which allow us to consider 
planetesimal disk formation in HD 69830 in the aftermath of planet 
formation. 

• We find that the model of lAlibert et all ( |2006|) in its basic 
form cannot coherently explain the observed eccentricities in the 
HD69830 system. Additional mass, either in the form of additional 
planets or increased masses of the observed planets, would need to 
be present in the system in order to excite greater eccentricities in 
the observed planets. 

• Some way of mitigating the effect of gas damping on the plan- 
ets or stimulating additional eccentricity excitation after the end of 
the gas damping phase would also seem to be required. 

• The eccentricities of the planets are relatively tmimportant in 
deciding the overall availability of an excited planetesimal popula- 
tion with pericentre(s) outside 1 AU. 

• We find that the eccentricity damping of planetesimals does 
not act as an insurmountable obstacle to the existence of an excited 
eccentric disk: We consistently find in all simulations that ~ 15% 
(equivalent to ~ 25Me) or more of the total solid disk material will 
remain in excited orbits having pericentres, ^ > 1 AU. It therefore 
seems probable that after the planet formation process had con- 
cluded, HD 69830 would have been left with a significant swarm 
of eccentric planetesimals. 

• We cannot yet definitively rule-out any of the IR emission 
models discussed in 91.21 but the consistent production and sur- 
vival of eccentric disks during and after the planet formation pro- 
cess in our n-body models does suggest the eccentric swarm idea is 
worthy of further investigation. 

• Gas damping of planetesimals does not significantly alter the 
number of available planetesimals, but it does serve to significantly 
decrease the proportion of these planetesimals which have high ec- 
centricities (e > 0.5). 
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• Gas damping of planetesimals works to "size-sort" the plan- 
etesimals, preferentially leaving the larger planetesimals (which 
contain a large proportion of the mass) occupying the higher eccen- 
tricity orbits whilst the low mass objects are efHciently circularised. 

• Further work is needed (and is in progress) to model the life- 
time and emission characteristics of extended non-aligned eccen- 
tric swarms of planetesimals to understand in more detail whether 
the reservoirs of planetesimals predicted in this paper could survive 
long enough to explain the observed IR emission. 
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